Analytic study of the urn model for separation of sand 
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We present an analytic study of the urn model for separation of sand recently introduced by 
Lipowski and Droz (Phys. Rev. E 65, 031307 (2002)). We solve analytically the master equation and 
the first-passage problem. The analytic results confirm the numerical results obtained by Lipowski 
and Droz. We find that the stationary probability distribution and the shortest one among the 
characteristic times are governed by the same free energy. We also analytically derive the form of 
the critical probability distribution on th e critical line, whic h supports their results obtained by 
numerically calculating Binder cumulants ( sond-mat /020 1 472| ) . 
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I. INTRODUCTION 

A granular system exhibits extremely rich phenomena, 
which has recently attracted extensive studies. One of 
such interesting phenomena is the spatial separation of 
shaken sand Sand in a box separated into two equal 
parts by a wall that allows the transfer of sand through 
its narrow slit prefers to aggregate more in one side under 
certain conditions. 

Eggers explained the emergence of symmetry breaking 
using a hydrodynamic approach [Q. The key idea is to 
introduce the effective temperature taking into account 
the inelastic collisions for granular material. 

Lipowski and Droz proposed a dynamic model to ex- 
plain the essence of the phenomena The model 
is a certain generalization of the Ehrenfest model [Q. 
Interestingly this model shows a spontaneous symme- 
try breaking in contrast to other generalizations of the 
Ehrenfest model. They derived the master equation and 
found in a numerical way the phase diagram that displays 
a rich structure like continuous and discontinuous transi- 
tions as well as a tricritical point. They also numerically 
solved the first-passage problem to find exponential or 
algebraic divergences. 

Thanks to its simplicity, the model allows analytic ap- 
proaches. In this paper, we present the results of this an- 
alytic study to the master equation and the first-passage 
problem addressed by Lipowski and Droz. These not only 
confirm their numerical results but also give us some in- 
sights in the nature of the discontinuous transition in the 
stationary probability distribution. We also analytically 
derive the form of the critical probability distribution on 
the critical line. 

The paper is organized as follows. In Sec. II we brie fly 
review the model and its master equation. In Sec. Ill we 



present the analytic solution of the master equation in 
the thermodynamic limit and the analytic expression of 
the stationary probability distribution. The form of the 
stationary probability distribution on the critical line is 
also derived. In Sec. IV we analytically solve the first- 
passage problem. Detailed analysis on the behavior of 



the characteristic times is given in Sec. Section VI is 
devoted to conclusions and discussions. 



II. MODEL AND ITS MASTER EQUATION 

The model introduced by Lipowski and Droz |^ is de- 
fined as follows. N particles are distributed between two 
urns, and the number of particles in each urn is denoted 
as M and N — M , respectively. At each time of up- 
dates one of the N particles is randomly chosen. Let n 
be a fraction of the total number of particles in the urn 
that the selected particles belongs to. With probability 
exp(— the selected particle moves to the other urn. 
T{n) represents the effective temperature of an urn with 
particles nN that measures the thermal fiuctuations of 
the urn. Lipowski and Droz chose the temperature as 
r(n) =To-h A(l-n). 

It is easy to derive the master equation for the prob- 
ability distribution p(M, t) that there are M particles in 
a given urn at time t S 



p{M,t+l) = 



+F{^)p{M+l,t) 



N 
1 - F 



where F{n) — nexp(— j^;^) measures the flux of parti- 
cles leaving the given urn. Here we introduced for con- 
venience the notations p(— 1, t) = p{N -|- 1, t) = 0. 

The difference in the occupancy of the urns can be 
represented by the particle excess ^ = ^ 5- The time 
evolution of the averaged particle excess e{t) — (e)t = 
Ea/Itt ~ l )Pi^^^ i) is governed by 



(2) 



where !F{e) = F{^ — e) — F{^ + e) measures the net flux of 
particles in the given urn. One conventionally takes the 
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FIG. 1: Phase diagram of the urn model The symmet- 
ric solution vanishes continuously on the solid line while the 
asymmetric one disappears discontinuously on the dotted line. 
The transition of the behavior of the stationary probability 
distribution is denoted by the dashed line. 



unit of time in such a way that there is a single update 
per a particle on average. Therefore we scale the time by 
A''. Expanding Eq. (||) with respect to and using the 
mean-field approximation in evaluating the average, we 
get 



FIG. 2: A(e) for A = 0.3, To = 0.2. The mapping for e > 
is analytically continuated. 



Scaling again the time by N, expanding Eq. (Q) with 
respect to and noting that the second term in the 
right-handed side can be combined into a total derivative 
with respect to e, we finally obtain the partial differential 
equation 



0. 



(5) 



(3) 



Note that the stationary solution of Eq. (||) is determined 
by zero points of .?-"(•). The stable stationary solutions are 
given by zero points of J^{-) with a negative slope, which 
we call as the stable fixed points while the unstable ones 
are given by those with a positive slope, which we call as 
the unstable fixed points. 

Detailed analysis on the existence of the stable station- 
ary solutions of Eq. (^) was done by Lipowski and Droz 
We here display their phase diagram in Fig. Q to 
make our paper as self-contained as possible. The stable 
symmetric solution (e — 0) exists in region I, III, and 
IV while the stable asymmetric solution (e > 0) exists in 
region II, III and IV. 



III. THE SOLUTION OF THE MASTER 
EQUATION 

We are mainly interested in investigating the proper- 
ties of the infinite system. Consider the thermodynamic 
limit N ^ oo with ^ = 5 + e being fixed. Represent- 
ing the probability distribution by e instead of M, and 
expanding Eq. (|l|) with respect to and keeping the 
terms up to the first order, we arrive at the expression 



p(e,t+l) = p{e,t) + 



1 

N 



(4) 



Note that is zero at a finite number of points. The 
solution of Eq. (|^) can be found in the intervals that do 
not include those points. At each interval, it would be 
convenient to introduce a new variable 



A(6) 



dx 
0^ 



(6) 



where eg is a certain point in the interval. Figure g shows 
a typical behavior of the mapping. We also displayed the 
map for |e| > i by analytic continuation. This is nec- 
essary since the solution of Eq. (||) is of wave nature. 
(See Eq. (H) below.) A increases as e approaches to the 
stable fixed points while it decreases as e approaches to 
the unstable fixed points. With the help of this parame- 
terization and setting R{X,t) — T{e)p{e,t), Eq. now 
takes the form 



-i?(A,t) + — i?(A,t) = 0. 



(7) 



Note that Eq. (0) is in fact a half part of the wave equa- 
tions so that its solution is written as R{X,t) = /(A — t) 
with /(•) being an arbitrary differentiable function. It 
represents a wave that moves to the direction of increas- 
ing A as time evolves, which means that the system moves 
to stable fixed points. The solution of the original partial 
differential equation (|^) now reads 



(8) 
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From the initial probability distribution po{e) = p{e, 0) 
we can determine the function /(•). So we get 



where et is given by the relation 

A(et) = A(e)-i. 



(9) 



(10) 



Here it should be understood that et is to be chosen in the 
same interval where e belongs to. Furthermore po(e) = 
for |e| > ^ is assumed since Ct in Eq. ^ can be larger (or 
smaller) than i (-5)- This happens because of the nature 
of the wave solution. Using the mapping from e to et, it 
is straightforward to show that the total probability is 
conserved: 



p{e,t)de 



po{et)det = 1 . 



(11) 



The shape of the probability distribution is distorted 
by a ratio !F{et)/^{(-) so that it accumulates at the near- 
est stable fixed points. In Eq. (^), the ratio approaches 
zero as time evolves unless e is on a stable fixed point. 
As a consequence, in the long time limit t — s- 00 the prob- 
ability distribution becomes a sum of delta peaks at the 
stable fixed points ti 



^(e.oo) ^^p^6{e- ei) , 



(12) 



where pi are the sum of the initial probabilities in two in- 
tervals adjacent to its point a. We would like to point out 
that the system is not ergodic and its dynamical phase 
space is decomposed into disconnected sectors. Each sec- 
tor is associated with a stable fixed point and is separated 
by the unstable fixed points. 

The fixed point condition !F{e) = is equivalent to 
Eq. (4) in Ref. where Lipowski and Droz analyzed 
in detail the condition and their results are summarized 
in the phase diagram (See Fig. |l|.). We would like to 
mention that in regions III and IV in Fig. |l| both the 
symmetric and the asymmetric solutions exist together. 
In fact, either solution can be realized by choosing an ap- 
propriate initial configuration. Lipowski and Droz distin- 
guished regions III and IV according to the different be- 
haviors of the stationary probability distribution in their 
numerical process of taking the limit N 00. In region 
III there appear two delta peaks for the asymmetric so- 
lutions while in region IV there appears only the central 
delta peak for the symmetric solution. It is contradic- 
tory to our result Eq. ( p^ ) where any delta peaks for the 
stable fixed points can appear depending on the initial 
configurations. 

To resolve this contradiction and understand the na- 
ture of the transition between regions III and IV, we take 
another limit in the master equation (|l|) , namely take the 
long time limit t 00 before we take the limit N ^ 00. 



This limit may not properly reflect the properties of the 
infinite system. Since the infinite system is not ergodic 
as we showed above, changing the order of taking lim- 
its iV —> cxd and t 00 may not yield the same result. 
Anyway it seems that in their simulations about the sta- 
tionary probability distribution, Lipowski and Droz took 
the limit t — > 00 for a finite system size N, and then 
extrapolating the results to iV ^ 00. 

Let's first take the long time limit of i ^ 00 in Eq. (|^). 
In this limit we may drop off the time dependence in the 
probability distribution, which now takes the form 



p{N) 
p(M) 



F{^)p{N~\) + {l-F(l))p{N) 
+F(^^)p{M + l) 



p(0) 



N 

for M - 1. 
1 



N -M^ 
2,1 



piM) 



F{-)pil)+{l~F{l))piO). 



(13) 



The first equation in Eq. ( |13[ ) allows us to rewrite p{N) 
in terms of p{N — 1), which in turn allows to rewrite 
p{N — 1) in terms of p{N — 2), and so on. Therefore we 
find 



p{M) = 



■p{M - 1) 



p(o)n 



F 



N 



(14) 



for M = N, . . . ,2,1. p{0) appears as an overall factor to 
normalize the probabilities so that we get 

N M 

p{0) = ' 



1 



En 

M=l i=l 



F 



(15) 



Now let's take the limit iV ^ 00. With ^ = ^ + e, and 
jj: = 5 + 2;, and scaling the probability distribution by N, 
the stationary probability distribution for large N now 
becomes 



Psie) 



,AfG(€) 



(16) 



■'-2 
where 

G(e) - l\ dx [logF(i - x) - logF(i + x)] . (17) 

In the limit N 00, the main contribution to the 
stationary probability distribution comes only from the 
maximum of G{-), and it becomes delta peaks. The max- 
imum of G(e) occurs when 



G'(e)=logi^( 



1 



logF( 



1 



0. 



(18) 
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Since G"(e) is the difference of logarithms of F(i — e) and 
+ e), both G{-) and JF(-) share the similar quahta- 
tive properties. For example, the maximum of both func- 
tions occurs at the same stable fixed points. Note that 
in region II only the two asymmetric solutions with pos- 
itive and negative particle excesses are stable, and have 
the same maximum while in region I only the symmetric 
solution is stable. Therefore the stationary probability 
distribution has the double peaks in region II and only 
the central peak in region I. In region III and IV both 
the symmetric and the asymmetric solutions are stable so 
that the maximum of G{e) should be determined by com- 
paring its values at the stable fixed points. The crossover 
of the maximum point occurs when both values coincide. 
This implies that the transition between the double peaks 
and the central single peak in the probability distribution 
is determined by the condition AG — G{ea) — G(0) — 0, 
where is the nonzero stable fixed point. This condition 
yields a line that separate two regions III and IV. 

It is very interesting to see that — G(-) resembles the 
free energy of the equilibrium systems, and the transi- 
tion between two phases is determined by the condition 
that the free energies of both phases are equal. Further- 
more a certain characteristic time behaves differently in 
two phases, as Lipowski and Droz numerically found 
We will show an analytic relation between them in next 
section. 



IV. CHARACTERISTIC TIMES 



is reached, we get the expression 



Ar(M) 



N-AI 

E 

1=1 j=i 



I M+j ' 
\ N , 



(21) 



Since t{^] 



by definition, we immediately get -t- 
1) = Ar(f + 1), which is given by Eq. (|2l|) with M = 
^4-1. By successively adding At(M), we get the general 



expression for t(M) for M — N, N — 1, 



N 
' 2 



1: 



r{M)= J2 



^ F(Jl) 
fc=". + l ^ \n) 



N-k i 

En 



F 



■ jV-Af-j + l - 
N 



• M +j ' 
. N , 



(22) 

We are mainly interested in the behavior of r(M) as 
N increases. Again we scale the characteristic times by 

iV. With # = 5 
1 

2 



e being fixed, and introducing = 
' X, the summations for large iV 
can be replaced by integrations so that the characteristic 
times for e > takes the form 



' N 



y- 



r(e) « N 



dx 



dye 



NH{x,y) 



(23) 



with H{x,y) — G{y) — G{x). The longest characteristic 
time t{N) corresponds to T(e — i). The shortest one 
t( Y + 1) corresponding to r(e = 0) is, in general, smaller 
than r(e) with positive e by an order of magnitude. It is 
necessary to deal with it separately. We get 



Lipowski and Droz defined an averaged first-passage 
time t{M) needed for a configuration with M particles 
in an urn (and N — M particles in the other urn) to 
reach the symmetric configuration (ill — ^) [H- They 
obtained the relations among the averaged characteristic 
times from the dynamical rules as 

r{M)^F{^)[T[M-l) + l] 

+ (19) 



Y + 1 • Here it is understood that 



forM = 7V,A^-l, 

the term associated with t{N + 1) does not appear. (In 
fact, its coefficient -F(O) vanishes.) By definition of the 
characteristic times, t(y) =0 and t{N — M) = t{M). 

Defining the difference of successive characteristic 
times as At(M) = t(M) - t(M - 1), Eq. (|l9|) can be 
rewritten as 



At(M) 



1 



l + F{^^L_E)AriM + l)\. (20) 



By applying this relation repeatedly until At{N) 



^ 2 



1) 



dye 



NH{0,y) 



(24) 



Since H{Q,y) ~ G{y) — G(0), both the shortest char- 
acteristic time t(-^ -|- 1) and the stationary probability 
distribution ps (e) have essentially the same functional de- 
pendence on G(-). Therefore it is not surprising that be- 
haviors of both quantities for large N are closely related. 
However it is not clear why they are. 



V. ANALYSIS OF THE CHARACTERISTIC 
TIMES 

We first consider the behavior of r( ^ + 1), the shortest 
one among the characteristic times. For large N, the 
main contribution to + 1) comes from the maximum 
point y,n of iJ(0, y) = G{y) — G(0), or G(y), which was 
dealt in Sec. IH. 



FTTT 



In region I and IV the maximum occurs at j/„i — 
corresponding to the symmetric configuration, while in 
region II and III it occurs at ym > corresponding to 
the asymmetric configurations. Therefore the maximum 
is zero in region I and IV while it is positive in region II 
and HI. 
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When j/„ 
to get 



0, we may expand H{0, y) around y = 



-2 

4 
~3 



1 - 
1 - 



A/2 



(To + A/2)2 
3(A/2)2 



(To + A/2)4 



(25) 



Therefore it yields + 1) ^ N^^ as long as the co- 
efficient of the first term in Eq. ( |25| ) is negative. This 
is the case in region I and IV. The coefficient vanishes 

when To — ~ y i which corresponds to the critical 

line found by Lipowski and Droz |^ . On this line, we get 



i?(0,y)«-^(l-^A)y4. 



32 
15 



(26) 



100 




FIG. 3: characteristic time t{N) as a function of N for 
To = 0.2 and A = 1.72, 1.744067, 1.77, 1.8, 2.0 (from the top). 



which yields r(-| + 1) ~ iV 



for A < |, andT(f + 1) 



N~B at A = I corresponding to the tricritical point. 

When Um > 0, the maximum is positive so that r( + 
1) ~ N~^e°'^ (with a being a positive constant), that 
is, it diverges exponentially. 

Now let's investigate the behavior of the longest char- 
acteristic time t{N) or r(e — 1). Again the main con- 
tribution comes from the maximum point of H{x,y) in 
region restricted by three lines y ~ x,x ^ 0,y = i. Note 
that y > X in the region. Interestingly the maximum 
point {xjn, y„i) is closely related with the fixed points of 
He). 

In region I, e = is only the fixed point so that x„i — 
ym = 0, and the maximum is zero. Expanding H{x,y) 
about this point, we get 



H{x,y) 



1 - 



A/2 



(To + A/2)2 



(j/'-.x^) + 



(27) 



The coefficient of the leading term in Eq. 
only if To > 



2 ' 



is negative 



that is, above the critical line. 



Changing the variables a;, y to the polar coordinates r, 9 
and scaling the radial coordinate r by \/]V, we arrive at 



t(N) 



drr dO exp 



-21 



A/2 



(To + A/2)2 
X cos(26')j . (28) 



Here c is a constant that gives the upper bound of the 
integration over r. The contribution from the neighbor- 
hood of = J yields a logarithmic divergence. Fig. ^ 
shows a typical behavior of characteristic time t{N) as a 
function of N for several values of A with Tq = 0.2. The 
first two uppermost lines stand for t{N) in region II, and 
the others represent that in region I. We conclude that 
in region I t{N) diverges logarithmically as N increases. 



As we see in Eq. (|2^) , the leading term vanishes on the 
critical line. On this line we need to expand more. So 



H{x,y) 



■i' 



(29) 



Again changing the variables x, y to the polar coordinates 
and scaling the radial coordinate appropriately (by Nt- 
or A^e), we conclude that t{N) diverges algebraically as 
A^2 on the critical line and as A^3 at the tricritical point. 

In regions II, III, and IV there appear many fixed 
points among which we can always find one with y,„ > 
Xra and G{]jm) > G{xm)- There, the maximum 
H{xm,yTn) is positive, and t{N) diverges exponentially 
as A^ increases. We would like to point out that the 
situation is different from that of r(-^ -I- 1). In fact, it 
corresponds to the case with Xm being fixed to zero. 

Finally let's consider the behavior of t{N) on the line 
separating two regions I and IV. As we approach this 
line from region IV, the nonzero fixed points merge to 
disappear at e = ei > 0. That is, we can find a positive ei 
such that J^i^i) — T' {ei) — 0. On this line the maximum 
point is given by Xm = Vm = ei so that the maximum 
H{xjm y-m) is zero. Expanding H{x, y) around this point, 
we get 



H[x,y) 



G"'{e,){{y-e,f~{x-e,f)+. 



(30) 



(Here we don't write down G""(ei) explicitly since it is 
not important as far as it is positive.) Note that the 
leading order is the third instead of the fourth as in case 
of the critical line. The reason is that G{() is not sym- 
metric about e — ei while it is symmetric about e = 0. 
Consequently t{N) on this line diverges algebraically as 
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VI. CONCLUSIONS 

We analytically investigate the urn model introduced 
by Lipowski and Droz g]. We exactly solve the mas- 
ter equation of the model in the thermodynamic limit 
and find how the probability distribution evolves. In 
the long time limit, the probability distribution becomes 
delta peaks only at the stable fixed points. In fact the 
ergodicity of the dynamics is broken so that the dynam- 
ical phase space is decomposed into disconnected sectors 
separated by the unstable fixed points. The strength of 
a delta peak is equal to the sum of initial probabilities in 
the disconnected sector it belongs to. 

We also solve exactly the stationary probability distri- 
bution where we take the long time limit before we take 
thermodynamic limit. Regardless of the initial probabil- 
ity distribution it shows double peaks or a single cen- 
tral peak depending on the parameters of the system. 
The final formula of the stationary probability distribu- 
tion resembles that of the equilibrium systems, where 
the transition from the doubles peaks to the single peak 
is determined by the condition that free energies of two 
phases become equal. 

Recently Lipowski and Droz numerically calculated 
Binder cumulants of the urn model to find that the criti- 



cal probability distribution has the form p{x) ^ e ^ on 
the critical line, and p{x) ~ e~^ on the tricritical point, 
where x is the rescaled order parameter proportional to 
the particle excess e. As we showed, G(e) ^ i?(0, e) and 
its behavior on the critical line (including the tricritical 
point) is given by Eq. (p6|). Therefore the critical proba- 
bility distribution has the form p(e) ~ exp[— — |A)e^] 
on the critical line, andp(e) ~ exp[— yUe^] on the tricrit- 
ical point. Our analytic result supports their numerical 
result. 

The first-passage problem is analytically solved. In- 
terestingly both the shortest characteristic time and the 
stationary probability distribution are governed by the 
same free energy. Therefore the behavior of the shortest 
characteristic time and the properties of the stationary 
probability distribution are closely related. The analytic 
results on the behavior of the characteristic times sup- 
port the numerical results of Lipowski and Droz |^ . 

Finally it would be very interesting to understand 
why both the stationary probability distribution and the 
shortest characteristic time are governed by the same free 
energy. It would also be of interest to extend our ana- 
lytic study to many-urn models and other types of urn 
models. 
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